{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Sampling and Hypothesis Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hypothesis Tests | Parametric tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parametric tests rely on a probability distribution of known form as a model for the null hypothesis. \n",
    "\n",
    "Here we will look at some commonly-encountered examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How surprising is my result? Calculating a p-value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many circumstances where we simply want to check whether an observation looks like it is compatible with the null hypothesis, $H_{0}$. \n",
    "\n",
    "Having decided on a significance level $\\alpha$ and whether the situation warrants a one-tailed or a two-tailed test, we can use the cdf of the null distribution to calculate a p-value for the observation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: probability of rolling a six\n",
    "\n",
    "Your arch-nemesis Blofeld always seems to win at ludo, and you have started to suspect him of using a loaded die.\n",
    "\n",
    "You observe the following outcomes from 100 rolls of his die:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "data <- c(6, 1, 5, 6, 2, 6, 4, 3, 4, 6, 1, 2, 5, 6, 6, 3, 6, 2, 6, 4, 6, 2,\n",
    "       5, 4, 2, 3, 3, 6, 6, 1, 2, 5, 6, 4, 6, 2, 1, 3, 6, 5, 4, 5, 6, 3,\n",
    "       6, 6, 1, 4, 6, 6, 6, 6, 6, 2, 3, 1, 6, 4, 3, 6, 2, 4, 6, 6, 6, 5,\n",
    "       6, 2, 1, 6, 6, 4, 3, 6, 5, 6, 6, 2, 6, 3, 6, 6, 1, 4, 6, 4, 2, 6,\n",
    "       6, 5, 2, 6, 6, 4, 3, 1, 6, 6, 5, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do you have enough evidence to confront him?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will work with the binomial distribution for the observed number of sixes\n",
    "\n",
    "# Write down the hypotheses\n",
    "# H0: p = 1/6\n",
    "# H1: p > 1/6\n",
    "\n",
    "# choose a significance level\n",
    "# alpha = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  [1]  TRUE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE\n",
      " [13] FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE\n",
      " [25] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE\n",
      " [37] FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE\n",
      " [49]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE\n",
      " [61] FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE\n",
      " [73] FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE\n",
      " [85]  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE\n",
      " [97]  TRUE  TRUE FALSE FALSE\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "43"
      ],
      "text/latex": [
       "43"
      ],
      "text/markdown": [
       "43"
      ],
      "text/plain": [
       "[1] 43"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "100"
      ],
      "text/latex": [
       "100"
      ],
      "text/markdown": [
       "100"
      ],
      "text/plain": [
       "[1] 100"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# code the data as 6=success and {0-5}=failure\n",
    "six <- data==6\n",
    "print(six)\n",
    "\n",
    "# how many sixes were observed?\n",
    "x <- sum(six)\n",
    "x\n",
    "\n",
    "# check number of trials\n",
    "n <- length(data)\n",
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "5.43908695860296e-10"
      ],
      "text/latex": [
       "5.43908695860296e-10"
      ],
      "text/markdown": [
       "5.43908695860296e-10"
      ],
      "text/plain": [
       "[1] 5.439087e-10"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now use H0 to find the p-value of the observed number of sixes\n",
    "pval <- 1 - pbinom(42,100,1/6)  # note this uses k=(observed value-1)\n",
    "pval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pval is less than alpha, so reject H0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: is the coin fair?\n",
    "\n",
    "Dr Vogel has challenged you to a game of ludo and you have agreed to flip a coin to decide who starts.\n",
    "\n",
    "You're not sure whether the coin she is using is fair or not.\n",
    "\n",
    "She flips it 50 times for you, with the following results: (1=heads, 0=tails)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data <- c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does the coin appear to be fair? If not, should you pick heads or tails?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Again, this is a binomial distribution for the observed number of heads, but this time the test is 2-tailed\n",
    "\n",
    "# Write down the hypotheses\n",
    "# H0: p = 1/2\n",
    "# H1: p != 1/2\n",
    "\n",
    "# choose a significance level\n",
    "# alpha = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "20"
      ],
      "text/latex": [
       "20"
      ],
      "text/markdown": [
       "20"
      ],
      "text/plain": [
       "[1] 20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "50"
      ],
      "text/latex": [
       "50"
      ],
      "text/markdown": [
       "50"
      ],
      "text/plain": [
       "[1] 50"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find the number of heads\n",
    "h <- sum(data)\n",
    "h\n",
    "\n",
    "# check number of trials\n",
    "n <- length(data)\n",
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "25"
      ],
      "text/latex": [
       "25"
      ],
      "text/markdown": [
       "25"
      ],
      "text/plain": [
       "[1] 25"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# now use H0 to find the p-value of the observed number of heads\n",
    "ex <- 50 * 0.5 # the expected value\n",
    "ex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0.2026388\n"
     ]
    }
   ],
   "source": [
    "x1 <- 20  # the lower tail\n",
    "p1 <- pbinom(x1,50,0.5)  \n",
    "pval <- 2 * p1 # double the p-value for a two-tailed test\n",
    "print(pval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pval is greater than alpha, so we accept H0: there is no evidence that the coin is biased, at the 5% level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Difference between two means: independent 2-sample t-test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the **t test** to assess whether two samples taken from normal distributions have significantly different means. \n",
    "\n",
    "The test statistic follows a Student's t-distribution, provided that the variances of the two groups are equal.\n",
    "\n",
    "Other variants of the t-test are applicable under different conditions.\n",
    "\n",
    "The test statistic is\n",
    "\n",
    "$$ t = \\frac{\\bar{X}_{1} - \\bar{X}_{2}}{s_p \\cdot \\sqrt{\\frac{1}{n_{1}} + \\frac{1}{n_{2}}}} $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ s_p = \\sqrt{\\frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2}} $$\n",
    "\n",
    "is an estimator of the pooled standard deviation.\n",
    "\n",
    "Under the null hypothesis of equal means, the statistic follows a Student's t-distribution with $(n_{1} + n_{2} - 2)$ degrees of freedom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: difference in birth weight\n",
    "\n",
    "The birth weights of babies (in kg) have been measured for a sample of mothers split into two categories: nonsmoking and heavy smoking.\n",
    "\n",
    "- The two categories are measured independently from each other. \n",
    "- Both come from normal distributions\n",
    "- The two groups are assumed to have the same unknown variance.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_nonsmoking <- c(3.99, 3.79, 3.60, 3.73, 3.21, 3.60, 4.08, 3.61, 3.83, 3.31, 4.13, 3.26, 3.54)\n",
    "data_heavysmoking <- c(3.18, 2.84, 2.90, 3.27, 3.85, 3.52, 3.23, 2.76, 3.60, 3.75, 3.59, 3.63, 2.38, 2.34, 2.44)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to know whether there is a significant difference in mean birth weight between the two categories.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write down the hypotheses\n",
    "# H0: there is no difference in mean birth weight between groups: d == 0\n",
    "# H1: there is a difference, d != 0\n",
    "\n",
    "# choose a significance level\n",
    "# alpha = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "'non-smoking: n = 13 , mean = 3.66769230769231 , SD = 0.297941656653162'"
      ],
      "text/latex": [
       "'non-smoking: n = 13 , mean = 3.66769230769231 , SD = 0.297941656653162'"
      ],
      "text/markdown": [
       "'non-smoking: n = 13 , mean = 3.66769230769231 , SD = 0.297941656653162'"
      ],
      "text/plain": [
       "[1] \"non-smoking: n = 13 , mean = 3.66769230769231 , SD = 0.297941656653162\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "'heavy smoking: n = 15 , mean = 3.152 , SD = 0.514006670329148'"
      ],
      "text/latex": [
       "'heavy smoking: n = 15 , mean = 3.152 , SD = 0.514006670329148'"
      ],
      "text/markdown": [
       "'heavy smoking: n = 15 , mean = 3.152 , SD = 0.514006670329148'"
      ],
      "text/plain": [
       "[1] \"heavy smoking: n = 15 , mean = 3.152 , SD = 0.514006670329148\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "n_ns <- length(data_nonsmoking)\n",
    "n_hs <- length(data_heavysmoking)\n",
    "\n",
    "mean_ns <- mean(data_nonsmoking)\n",
    "mean_hs <- mean(data_heavysmoking)\n",
    "\n",
    "s_ns <- sd(data_nonsmoking)\n",
    "s_hs <- sd(data_heavysmoking)\n",
    "\n",
    "paste(\"non-smoking: n =\", n_ns, \", mean =\", mean_ns, \", SD =\",s_ns)\n",
    "paste(\"heavy smoking: n =\", n_hs, \", mean =\", mean_hs, \", SD =\",s_hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.515692307692308"
      ],
      "text/latex": [
       "0.515692307692308"
      ],
      "text/markdown": [
       "0.515692307692308"
      ],
      "text/plain": [
       "[1] 0.5156923"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# difference between the two sample means:\n",
    "d_obs <- mean_ns - mean_hs\n",
    "d_obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.428057812829366"
      ],
      "text/latex": [
       "0.428057812829366"
      ],
      "text/markdown": [
       "0.428057812829366"
      ],
      "text/plain": [
       "[1] 0.4280578"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the pooled standard deviation\n",
    "sp <- sqrt(((n_ns - 1)*s_ns^2 + (n_hs - 1)*s_hs^2)/(n_ns + n_hs - 2))\n",
    "sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "3.17926343494134"
      ],
      "text/latex": [
       "3.17926343494134"
      ],
      "text/markdown": [
       "3.17926343494134"
      ],
      "text/plain": [
       "[1] 3.179263"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the test statistic\n",
    "t_obs <- d_obs/(sp * sqrt(1/n_ns + 1/n_hs))\n",
    "t_obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "26"
      ],
      "text/latex": [
       "26"
      ],
      "text/markdown": [
       "26"
      ],
      "text/plain": [
       "[1] 26"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# degrees of freedom is given by n1 + n2 - 2\n",
    "df <- n_ns + n_hs - 2\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "2.05552943864287"
      ],
      "text/latex": [
       "2.05552943864287"
      ],
      "text/markdown": [
       "2.05552943864287"
      ],
      "text/plain": [
       "[1] 2.055529"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find the critical value\n",
    "t95 <- qt(1-0.05/2,df) # critical value for 95% of probability mass\n",
    "t95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# t_obs lies outside the 95% confidence interval [-t95,t95], so we reject H0"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
